function [mx, my, mx_list, my_list] = mxmy(bx, by, t, x, z0, zc, varargin)
%M_XY Summary of this function goes here
%   Detailed explanation goes here
    p=inputParser;
    pTest=@(x) x >= 1;
    nTest=@(x) x >= 0;
    p.addParameter('pmax', 5, pTest);
    p.addParameter('nmax', 10, nTest);
    p.parse(varargin{:});
    
    mx_list = zeros( p.Results.pmax+1, length(t) );
    my_list = zeros( p.Results.pmax+1, length(t) );
    for p_index=0:p.Results.pmax
        [mx_p, my_p] = arrayfun(@(t) m_p(bx, by, t, x, z0, zc, p_index, p.Results.nmax), t);
        mx_list(p_index+1, :) = mx_p;
        my_list(p_index+1, :) = my_p;
    end
    mx = sum(mx_list, 1);
    my = sum(my_list, 1);
end

